#!/usr/bin/perl
# 2010/05/12 Listファイルから配列を取得

use strict;
use warnings;
use IO::File;
use Getopt::Std;
use Carp qw(croak);
$| = 0;

#---------------------------------------
# ID fileの設定は必須
# -s　オプションなしだと、coding領域をそのまま取得する
# -s オプションで指定された場所からStart部位までを取得する
#---------------------------------------
if((@ARGV != 2  and @ARGV != 4) or $ARGV[0] eq '--help'){
	die "Usage : $0 -i [ID list file] -s(if get outside) [length from translation start site which you want to get]\n";
}
my %opt;
getopts("i:s:", \%opt);
my $ListFile = $opt{i};
my $GETLength = 0;
if(exists $opt{s}){
	$GETLength = $opt{s};
}

### 配列を準備する
my %SQ;
my $Dir = '/home/hideki/Data/S.Pombe_seq/sanger_latest-20090706/';
for(my $i = 1; $i <= 3; $i++){
	my $chr = 'chr' . $i;
	my $FILE = $Dir . $chr . '.fa';
	my $fh = IO::File->new($FILE) or die "cannot open $FILE: $!";
	$fh->getline();
	while($_ = $fh->getline()){
		s/\r?\n//;
		$SQ{$chr} .= $_;
	}
	$fh->close();
}

### すべてのCDSの場所
my $ALLCDS = '/home/hideki/Data/S.Pombe_seq/CDS_2007.bed';
#chr1    7618    9274    SPAC212.09c     0       +       7618    9274    0       1       1656,   0,
#chr1    12157   12994   SPAC212.08c     0       +       12157   12994   0       1       837,    0,
#chr1    13664   14555   SPAC212.07c     0       +       13664   14555   0       1       891,    0,
#chr1    18041   18557   SPAC212.06c     0       +       18041   18557   0       2       265,209,        0,307,
my %data;
{
	my $fh = IO::File->new($ALLCDS) or die "cannot open $ALLCDS: $!";
	while($_ = $fh->getline()){
		s/\r?\n//;
		my ($chr, $start, $stop, $name, $value, $direction, @other) = split /\t/;
		$data{$name} = [$chr, $start, $stop, $direction];
	}
	$fh->close();
}

### リストを読み込んで配列を取得する
{
	my $missinNumber = 0;
	my $fh = IO::File->new($ListFile) or die "cannot open $ListFile: $!";
	while($_ = $fh->getline()){
		s/\r?\n//;
		my $id = $_;
		unless(exists $data{$id}){
			warn ("$id was not found\n");
			$missinNumber++;
			next;
		}
		my $chr = $data{$id}->[0];
		my $start = $data{$id}->[1];
		my $end = $data{$id}->[2];
		my $direction = $data{$id}->[3];
		
		if(exists $opt{s}){
			### 外側を取得する場合　$GETLength　bp 分だけコード領域の前を取得
			if($direction eq '+'){
				my $sequence = lc(substr($SQ{$chr}, ($start - 1 - $GETLength), $GETLength));
				print ">$id\n";
				print "$sequence\n";
			}elsif($direction eq '-'){
				my $sequence = lc(substr($SQ{$chr}, ($end - 1), $GETLength));
				(my $reverse = (reverse $sequence)) =~ tr/atgc/tacg/;
				print ">$id\n";
				print "$reverse\n";
			}
			
		}else{
			###　コード領域を取得する場合
			my $length = $end - $start + 1;
			my $sequence = lc(substr($SQ{$chr}, ($start - 1), $length));
			
			print ">$id\n";
			if($direction eq '+'){
				print "$sequence\n";
			}elsif($direction eq '-'){
				(my $reverse = (reverse $sequence)) =~ tr/atgc/tacg/;
				print "$reverse\n";
			}			
		}
	}
	$fh->close();
	if($missinNumber != 0){
		warn ("total $missinNumber couldn't obtain the sequence\n");
	}
}

